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ABSTRACT 

We perform numerical simulations of nonlinear MHD waves in a gravitationally stratified molecular 
cloud that is bounded by a hot and tenuous external medium. We study the relation between the 
strength of the turbulence and various global properties of a molecular cloud, within a 1.5-dimensional 
approximation. Under the influence of a driving source of Alfvenic disturbances, the cloud is lifted up 
by the pressure of MHD waves and reaches a steady-state characterized by oscillations about a new 
time-averaged equilibrium state. The nonlinear effect results in the generation of longitudinal motions 
and many shock waves; however, the wave kinetic energy remains predominantly in transverse, rather 
than longitudinal, motions. There is an approximate equipartition of energy between the transverse 
velocity and fluctuating magnetic field (as predicted by small-amplitude theory) in the region of the 
stratified cloud which contains most of the mass; however, this relation breaks down in the outer regions, 
particularly near the cloud surface, where the motions have a standing-wave character. This means 
that the Chandrasekhar-Fermi formula applied to molecular clouds must be significantly modified in 
such regions. Models of an ensemble of clouds show that, for various strengths of the input energy, the 
velocity dispersion in the cloud a oc Z - 5 , where Z is a characteristic size of the cloud. Furthermore, a 
is always comparable to the mean Alfven velocity of the cloud, consistent with observational results. 

Subject headings: ISM: clouds - ISM: magnetic fields - MHD - methods: numerical - turbulence - waves 



1. INTRODUCTION 

Interstellar molecular clouds, the sites of current star 
formation in our Galaxy, have long been known to yield 
| supersonic line widths of molecular spectral lines (e.g., 
see Zuckerman & Palmer 1974). Objects classified as 
molecular clouds span a large range of mean radii (R ~ 
1 pc- 100 pc), masses (M ~ 10 2 M - 10 6 M ) and mean 
number density (n ~ 10 1 cm" 3 — 10 3 cm -3 ). In fact, these 
quantities are correlated with the one-dimensional veloc- 
ity dispersion a through the well known line-width-size- 
density relations (e.g., Solomon et al. 1987): 

a = 0.72 (i?/pc) 5 kms -1 , (1) 
n = 2.3x 10 3 (#/pc) _1 cm -3 . (2) 

Thus, the velocity dispersion is typically supersonic since 
the sound speed c s is only w 0.2 km s -1 for the typical 
molecular cloud temperature T « 10 K (e.g., Goldsmith 
& Langer 1978). 

The largest clouds, of mass M J> 10 M , often classi- 
fied as Giant Molecular Clouds (CMC's) are in fact com- 
plexes of smaller clouds, since the volume averaged den- 
sity may be lower than the excitation density of CO (e.g., 
Blitz & Williams 1999), and also lower than allowed from 
thermal stability arguments (Falgarone & Puget 1986). 
Hence, the basic building blocks of interstellar molecu- 
lar clouds, which contain most of mass of molecular ma- 
terial, are the dark (or dwarf) molecular clouds, which 
have R ~ 1 pc - 10 pc, M ~ 10 2 M - 10 4 M and 
n ~ 10 2 cm -3 — 10 3 cm -3 . These clouds have velocity dis- 
persion (7^1 — 2 km s" 1 , and represent the class of objects 
that we are interested in modeling in this study. We also 
note that observed smaller scale (i? ~ 0.1 pc) dense cores 
are a separate object class which are embedded within dark 



clouds and collectively contain only a small fraction of the 
total cloud mass. 

While even dark clouds have masses that significantly 
exceed the thermal Jeans (1928) mass 



/ T \ 3/2 / n \-i/2 

= 17 (iok) (ijw) < 3 > 

(where we have used p — mn, and m = 2.33mn, in which 
«ih is the mass of a hydrogen atom), the line- width-size- 
density relations do imply that molecular clouds are indi- 
vidually in an approximate virial balance between turbu- 
lent and gravitational energies. In this paper, we equate 
the presence of nonthermal line widths with the presence of 
a random superposition of nonlinear (presumably hydro- 
magnetic) waves, and refer loosely to the latter as "turbu- 
lence". Within each cloud, the turbulence is expected to 
collectively exert a force (e.g., Chandrasekhar 1951) which 
resists the inward pull of gravity. 

The origin and persistence of turbulent motions in 
molecular clouds remain an active area of investigation. It 
was long understood (see Mestel 1965; Goldreich & Kwan 
1974) that supersonic hydrodynamic motions would de- 
cay rapidly through shocks, thereby creating a mystery of 
why the turbulence was commonly observed throughout 
the lifetime of molecular clouds. For GMC's, the lifetime 
is estimated to be a few x 10 7 yr (e.g., Blitz & Williams 
1999), which is at least a few times longer than the crossing 
time t c = 2R/a of the complexes. The smaller dark molec- 
ular clouds may have even longer lifetimes ~ 10 s yr (Shu, 
Adams. & Lizano 1987). Since slow and fast mode MHD 
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waves are also compressive and can be highly dissipative, it 
was suggested by Arons & Max (1975) that the transverse 
Alfven mode might be a long-lived component, thereby 
preventing a rapid overall collapse of the clouds and ex- 
plaining the observed persistence of the turbulence over 
their inferred lifetime. Mouschovias (1975) made the re- 
lated suggestion that the long-lived component may be due 
to standing waves, i.e., normal mode oscillations, of mag- 
netized clouds. Such global magneto-gravitational motions 
cannot be studied with a simple plane wave analysis in an 
infinite uniform medium. 

The magnetohydrodynamic (MHD) wave picture has 
been strengthened by the detection of large-scale magnetic 
fields within molecular clouds, through maps of polar- 
ized absorption and emission (e.g., Vrba, Strom, & Strom 
1976; Goodman et al. 1990; Schleuning 1998; Matthews 
& Wilson 2002), and Zeeman measurements of the line- 
of-sight magnetic field strengths (Crutcher 1999 and refer- 
ences therein). The latter imply that the magnetic energy 
(like the turbulent energy) is comparable to the gravita- 
tional potential energy; equivalently, the mass-to-flux ratio 
is close a critical value 



= Cl G- 1 ' 2 . 



(4) 



In the above equation, the constant C\ has been calcu- 
lated to be in the range 0.13-0.17 based on detailed two- 
dimensional equilibrium states (Mouschovias & Spitzer 
1976; Tomisaka, Ikeuchi, & Nakamura 1988) and is equal 
to l/(27r) in the case of a flattened one-dimensional layer 
(Nakano & Nakamura 1978). Collectively, the magnetic 
field measurements allow the possibility that clouds are 
supported lateral to the large-scale field B by its associ- 
ated Lorentz force, while these magnetic field lines act as 
a carrier of Alfvenic disturbances with SB ~ B, explain- 
ing the observed spectral line widths and preventing the 
clouds from assuming a very flattened configuration, 

The dynamical effect of propagating MHD waves using 
analytic or semi-analytic means has been studied by sev- 
eral authors. Dewar (1970) developed a formalism for cal- 
culating the effect of small amplitude hydromagnetic waves 
on a slowly varying background medium, using the WKB 
approach. Assuming ideal MHD and no dissipation of the 
waves, this leads to a steady-state relation between wave 
pressure P w and gas density p of the form P w oc p 1 / 2 (Mc- 
Kee & Zweibel 1995). Simplified calculations of the effect 
of small amplitude MHD waves on a molecular cloud have 
been performed by Fatuzzo & Adams (1993) and Martin, 
Heyvaerts, & Priest (1997). Both satisfy the above scaling 
of P w in the ideal MHD limit. In particular, the WKB 
model of Martin et al. (1997) yields a steady-state density 
structure of an infinite one-dimensional cloud supported 
by short- wavelength Alfven waves in the ideal MHD limit, 
and also when accounting for damping of linear waves by 
ion-neutral friction. In addition to ion-neutral friction, 
which damps even linear waves (Kulsrud & Pearce 1969), 
there are several nonlinear effects which will work to en- 
hance dissipation. The second order effect of a gradient 
in the magnetic pressure V SB 2 /8ir will in general lead 
to steepening of the waves followed by dissipation (Co- 
hen & Kulsrud 1974). Zweibel & Josafatsson (1983) state 
the form of this dissipation rate, which can dominate the 



process of ion-neutral friction for nonlinear and/or long 
wavelength modes. There are also other known nonlin- 
ear avenues for wave dissipation, such as the conversion of 
a parallel-propagating Alfven wave into an acoustic wave 
and another Alfven wave traveling in the opposite direc- 
tion (Sagdeev & Galeev 1969). 

Recently, several studies have resulted in a solution 
of the full set of nonlinear equations of ideal MHD us- 
ing finite difference approximations. Gammie & Ostriker 
(1996) performed a one-dimensional numerical simulation 
of MHD turbulence in a periodic domain. This has been 
followed up by several multi-dimensional simulations, also 
in a periodic domain (e.g., Stone, Ostriker, & Gammie 
1998; Ostriker, Gammie & Stone 1999; Mac Low et al. 
1998; Mac Low 1999; Padoan & Nordlund 1999; Ostriker, 
Stone, & Gammie 2001). These models impose prescribed 
velocity fluctuations, at the initial time and sometimes also 
throughout the computed time, and investigate the dissi- 
pation rate of turbulence, as well as various properties 
of turbulent fluctuations. One of the important results of 
these papers is that the decay time of the MHD turbulence 
is comparable to the crossing time over the driving scale 
of the turbulence. The nonlinear coupling of the Alfven 
to the fast and slow mode MHD waves is considered to be 
a significant source of the dissipation which reduces the 
lifetime of MHD turbulence relative to the ideal of pure 
Alfvenic turbulence. We note that periodic models repre- 
sent only a local region of a much larger molecular cloud, 
and maintain a fixed mean density. They cannot study 
global effects associated with the density stratification in 
a cloud or the existence of a cloud boundary. 

In this paper, we perform a different type of numerical 
simulation of MHD turbulence. We concentrate on one 
self-gravitating cloud and study the effect of the turbu- 
lence on the mechanical structure of the cloud. This cor- 
responds to an extension of the model studied by Martin et 
al. (1997) into a fully nonlinear counterpart. As an initial 
condition, we use a hydrostatic equilibrium between ther- 
mal pressure and self-gravity in a cloud that is bounded by 
an external high temperature medium. We input turbu- 
lent energy into the system continuously and sec how the 
mechanical equilibrium changes. Similar numerical simu- 
lations have been performed to study the propagation of 
Alfven waves in the solar chromosphere and corona (e.g., 
Hollweg, Jackson & Galloway 1982; Mariska & Hollweg 
1985; Hollweg 1992; Kudoh & Shibata 1999; Saito, Ku- 
doh, & Shibata 2001). We extend this class of model to a 
self-gravitating cloud, and study the relation between the 
strength of the turbulence and various global properties of 
a molecular cloud. This is the first of a series of papers on 
global models of MHD wave support in molecular clouds. 
In this paper, we develop high-resolution one-dimensional 
models under the assumption of ideal MHD. 

The paper is organized as follows. The numerical model 
we used for the simulation is summarized in § 2. The re- 
sults of the simulation are described in § 3. We add some 
discussion of the results in § 4 and summarize the paper 
in § 5. 



2. THE NUMERICAL MODEL 
2.1. Schematic View of our Model 
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Figure 1 shows a schematic picture of our model. We 
consider a molecular cloud that is threaded by a large- 
scale magnetic field, and concentrate on a local region of 
the molecular cloud enclosed by the rectangle in the fig- 
ure. We assume a driving force near the midplane of the 
cloud and follow the dynamical evolution of the vertical 
structure of the cloud. 
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without loss of generality, where v x and B x are the x- 
componcnts of the velocity and magnetic field, respec- 
tively. 

Therefore, the basic MHD equations we use in this paper 
are as follows: mass conservation, 



dp dp 



dv z 



the z-componcnt of the momentum equation, 



dv 



dv z 



1 dP 



1 B 9By +a ■ 
B„ — 1- g z , 



p dz Airp v dz 
the y-component of the momentum equation 



(8) 



(9) 



dv v 



+ v z 



8V V 



dt ' " z dz 4irp z dz 
the equation of energy, 



(10) 



Fig. 1. Schematic picture of our model. A molecular cloud 
that is threaded by a large-scale magnetic field is considered. 
A local region of the molecular cloud enclosed by the rectangle 
in the figure can represent the region we are modeling in this 
paper. We input a driving force near the midplane of the cloud 
and follow the dynamical evolution of the vertical structure of 
the cloud. 

2.2. Assumptions 

For simplicity, we assume ideal MHD, 1.5-dimensions, 
and isothermality for each Lagrangian fluid element. The 
1.5-dimensional approximation means that physical quan- 
tities depend on only one coordinate, but we evolve 
nonzero components of vectors in one additional direction. 
Isothermality for each Lagrangian fluid element means 
that the temperature does not change in time for each fluid 
element as it moves through Eulcrian space, which is dif- 
ferent from the assumption of a uniform time-independent 
temperature throughout the region. These assumptions 
appear in § 2.3 more concretely. 

2.3. Basic Equations 

We use local Cartesian coordinates (x,y, z) on the 
molecular cloud, where we set z to be the direction of 
the large-scale magnetic field. According to the symmetry 
of the 1.5-dimensional approximation, we set 



d_ _ d_ 

dx dy 



(5) 



The above symmetry and the divergence-free condition on 
the magnetic field imply 



B z = constant, 



(6) 



where B z is the z-component of the magnetic field that 
threads the molecular cloud. Moreover, from the assump- 
tion of linear polarization of the waves, we can set 



dT dT , ^^dv z 
the y-component of the induction equation, 



^By _ 

dt 

the equation of state, 



d 

— (-v z B y + v y B z ); 



where 



p kT 2 



Cs = 




(11) 



(12) 



(13) 



(14) 



is the isothermal sound speed; and the Poisson equation, 

(15) 



d 9z . n 



In these equations, t is the time, G is the gravitational con- 
stant, k is Boltzmann's constant, m is the mean molecular 
mass, 7 is the specific heat ratio, p is the density, P is the 
pressure, T is the temperature, g z is the ^-component of 
the gravitational field, v z is the z-componcnt of velocity, v y 
is the y-component of velocity, and B y is the y-componcnt 
of the magnetic field. 

In equation (11), we assume 7 = 1, so that the energy 
equation becomes 



dT ^ dT _ Q 

dt z dz 



(16) 



Br. = 



(7) 



which quantifies the assumption of isothermality for each 
Lagrangian fluid element. 
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2.4. Initial Conditions 

As an initial condition, we assume hydrostatic equilib- 
rium of a self-gravitating onc-dimcnsional cloud. The hy- 
drostatic equilibrium is calculated from the equations 



1 dP 

p dz 



= 9z 



dg z 
dz 



and 



= -4nGp, 



kT 



P = p- 

m 

subject to the boundary conditions 

g z (z = 0) = 0, 
p(z = 0) = po, 



and 



P(z = 0) = Po 



kT 



m 



(17) 
(18) 

(19) 

(20) 
(21) 

(22) 



where po and To are the initial density and initial temper- 
ature at z = 0, respectively. 

In order to solve the above equations, we need to assume 
an initial temperature distribution. If the temperature is 
uniform throughout the region, we have the following an- 
alytic solution ps found by Spitzer (1942); 



Ps(z) = Posech 2 (z/H ), 



where 



CsO 



is the scale height, and 



kT 
m 



(23) 
(24) 

(25) 



However, an isothermal molecular cloud is usually sur- 
rounded by warm or hot material, such as neutral hydro- 
gen or ionized gas. Therefore, we assume the initial tem- 
perature distribution to be 



T(z) = T + \{T C - To) 



1 + tanh f 1^ 



(26) 



where we take T c = 100Xo> z c = 3H , and Zd = 0.2H 
throughout the paper. This distribution shows that the 
temperature is uniform and equal to To in the region of 
< z < z c = 3H and smoothly increases to another uni- 
form value T c = 100T at z ~ z c = 3H - By using this 
temperature distribution, we can solve the ordinary dif- 
ferential equations (17)-(19) numerically. The numerical 
solution of these equations shows that the initial distribu- 
tion of density is almost the same as Spitzer's solution in 
the region of < z < z c (see Fig. 2). 

We also assume the following initial conditions; 

v z (z) =v y {z) =0, (27) 

B y (z) = 0, (28) 

B z (z) = B , (29) 

where B is a constant. According to equation (6), B z is 
spatially uniform and independent of time throughout the 
calculations. 



2.5. Driving Force 

We introduce a perturbation into the initially hydro- 
static cloud by adding a driving force, F(z,t), into the 
y-component of the momentum equation (10) as follows: 



pi ^ + Vz ^ ) ~^ Bz ^z- 



F(z,t), (30) 



where 



F(z,t) = 



and 



p a d(T<k) sin ( 27 "V) cx p[-(^) 2 ] 

pa d sin(27w i) exp[-(^) 2 ] 




to — Hq/ c s o. 



(t < m ) 

(10t < t < 40t ) 
{t > 40t ), 
(31) 

(32) 



Furthermore, is the amplitude of the induced accelera- 
tion, i>q is the frequency of the driving force, and z a repre- 
sents the region in which we input the driving force. The 
equations show that we input the sinusoidal driving force 
near the midplane of the cloud, and increase the maxi- 
mum driving force linearly with time until t — 10to, and 
maintain it to be constant during 10£o < t < 40£o- After 
t = 40£o, we terminate the driving force. 

2.6. Boundary Conditions 

We used a mirror symmetric boundary condition at 
z = and a free boundary at z — z out , the outer bound- 
ary of the calculation. In order to remove the reflection 
of waves at the outer boundary, we set z out to be a large 
value, i.e., 



z out ~ 2.6 x 10 4 ff , 



(33) 



and use nonuniform grid spacing for large z (see § 2.8). 

In addition to this, we set the gravitational field to be 
zero at large distances by introducing an artificial acceler- 
ation (g_) into equations (9) and (17), i.e., 



9-{z) = ~ 



9z{z) 



1 + tanh 



V 5#o 



(34) 



where we take z g = 1007? m this paper. This force is zero 
until z ~ z g = 100-f/o, but it becomes essentially equal to 
g z with negative sign for z > z g = 100i?o and compensates 
the gravity there. 

The above boundary conditions show that we have ef- 
fectively two outer boundaries. The first one is z ou t and 
the second one is z g , which is the boundary for the grav- 
itational field. In order to remove the reflection of waves 
at the outer boundary, it is useful to have a large z ut- 
However, ct v6ry large z ou t introduces numerical problems 
because the density decreases exponentially at large z, ac- 
cording to the stratification of a self-gravitating cloud. A 
very low density leads to a very large Alfven speed, making 
simulations inefficient by forcing very small time steps for 
an explicit calculation. By setting the net gravitational 
field to be zero at large distances, we can avoid an ex- 
tremely low density for large z and therefore take a large 
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z ou t in order to remove the reflection of Alfven waves at 
the outer boundary. We also note that for any slab of fi- 
nite extent, as opposed to the infinite slab implied by our 
one-dimensional model, the net gravitational field should 
indeed decrease at large distances, though not in the spe- 
cific manner prescribed here. 

Due to the added artificial acceleration, our effective 
numerical boundary is located at z ~ z g = 100i? , and 
we cannot trust the results beyond this region. However, 
as we show later, all dynamical events we are interested 
in occur within z < 50 Ho, where most of the mass and 
energy are concentrated. Hence, the effect of the artifi- 
cial acceleration is negligible for the main results in this 
simulation. 

2.7. Numerical Parameters 

A natural set of fundamental units for this problem are 
c s o, Ho, and po- These yield a time unit to — H$/ c s $. The 
initial magnetic field strength introduces one dimension- 
less free parameter, i.e. 



_ 8ttP 8np c 2 s0 
Po = —rjr ~ 



Bl 



(35) 



which is the initial ratio of gas to magnetic pressure at 
z = 0. 

In this cloud, /3q is related to the mass-to-flux ratio. 
For Spitzer's self-gravitating cloud, the mass-to-flux ratio 
normalized to the critical value is 



Bo' 



where 



S s = 



p s dz = 2p H 



(36) 



(37) 



is the column density of the Spitzer's self-gravitating 
cloud. Therefore, 



Po = Vs- 



(38) 



The column density of the cloud we used in this paper is 
almost equal to that of Spitzer's cloud. If we define the 
column density of the cloud, E, as the integral of density 
within — z c < z < z c , then 



p(t = 0)dz ~ 0.988 E s . 



(39) 



This means that we can use the value of /xs an excellent 
approximation to the dimensionless mass-to-flux ratio of 
the model cloud. 

Observations show that the dimensionless mass-to-flux 
ratios of molecular clouds are close to unity over a range 
of scales (Crutcher 1999; Shu et al. 1999). Hence, we take 
Po = 1 in the models presented in this paper. 

The driving force introduces three more important free 
parameters: a<j = ad{H / c 2 ) , the dimensionless ampli- 
tude of the acceleration due to driving, t> = the 
dimensionless frequency of driving, and z a = z a /H , the 
dimensionless scale of the driving region. For simplicity, 
we take vo = 1 and z a = 0.1 throughout this paper, and 
adjust the strength of driving by varying ad from 10 to 50. 

Dimensional values of all quantities can be found 
through a choice of T and po, along with the values of the 



dimensionless free parameters. For example, if To = 10 K 
and no = po/ m = f0 4 cm~ 3 , then c s o = 0.2 km s _1 , 
Ho = 0.05 pc, N s = T, s /m = 3xl0 21 cm" 2 , t = 2.5 xlO 5 
yr, and B = 20 if (3 = 1- 



2.8. Numerical Technique 

In order to solve the equations numerically, we use the 
CIP method (Yabc & Aoki 1991) for equations (8)-(10) 
and (16), and the MOCCT method (Stone & Norman 
1992) for equation (12). The combination of the CIP and 
MOCCT methods is summarized in Kudoh, Matsumoto 
& Shibata (1999). The CIP method is a useful method to 
solve an advection equation such as equation (16) accu- 
rately, and is also applicable to advection terms of equa- 
tions (8)-(10). 

Due to the mirror-symmetric boundary condition at 
z = 0, the Poisson equation (15) can be simply integrated 
from the midplane of the cloud; 



9z(z) 



AirG I p{z)dz. 
'o 



(40) 



We solve this equation by numerical integration. 

In this simulation, we actually use variants of the orig- 
inal CIP method. We use a conservative-CIP method for 
the mass conservation equation (8) , which was recently de- 
veloped by Xiao et al. (2002). This scheme assures exact 
conservation of mass and less numerical oscillations. The 
original CIP method does not assure the exact conserva- 
tion of mass, so that a systematic deviation during each 
time step can cause a problem for long time integrations; 
this is especially a problem if the Poisson equation, which 
utilizes the mass distribution, is being solved simultane- 
ously. Hence, the conservative-CIP method significantly 
improves the accuracy of our solution. We also use the 
monotonic-CIP method (Xiao, Yabe & Ito 1996) for the 
advection terms of equations (9), (10) and (16), and use 
the CCUP method (Yabe & Wang 1991) for the calcu- 
lation of gas pressure, in order to get more numerically 
stable results. The recent developments of CIP related 
schemes are summarized in Yabe, Xiao & Utsumi (2001). 

We used a uniform grid size, Azi = 0.02/Toj from z = 
to z = 80H , where Azi is the grid size of the i-th grid 
point. For 80-f/o < z < \2{)Ho, we gradually increase the 
grid size as Az.; + i = Min (1.05Azi, 0.2/fo), where Min is 
a function that picks the minimum from the two values; 
hence, the maximum grid size in this region is 0.2/fo- For 
120^0 < z < z out i where the gravitational field is compen- 
sated by the artificial force, we further increase the grid 
size according to Azi+\ — 1.05Azi. In this paper, we used 
a total of 4310 grid points, most of which arc concentrated 
in the region of the uniform grid, i.e., z < 80H - 



3. RESULTS 
3.1. Typical Result 
In this subsection, we show the results of a<j = 30 as the 



typical case. 
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3.1.1. Density Structure in the Cloud 

Figure 2 shows the density and temperature as a func- 
tion of z. The dashed lines show the initial distributions 
and the solid lines show the distribution at t = 30£o- 
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Fig. 2. Density (upper panel) and temperature (lower panel) 
as a function of z. The dashed lines show the initial distribu- 
tions and the solid lines show the distribution at t = 30io- The 
dotted line in the density plot shows the distribution of density 
in Spitzer's infinite self-gravitating equilibrium with a uniform 
temperature To. 

The dotted line in the density shows the distribution 
of Spitzer's self-gravitating cloud with uniform temper- 
ature of To. The initial density deviates from Spitzer's 
solution around z = z c = 3H where the initial tempera- 
ture increases up to 100Tb. The density decreases rapidly 
around z = z c = 3H because of the pressure balance 
between low temperature cloud and high temperature ex- 
ternal medium. For z > z c — 3H , the scale height of the 
density is very large and the density decreases gradually. 

The snapshot of density at t = 30io shows that the 
density has a complicated structure including many shock 
fronts. This is caused by the driving force in equation 
(30). The driving force generates nonlinear Alfven waves 
in the cloud which produce a magnetic pressure gradient. 
The magnetic pressure gradient and thermal pressure gra- 
dient usually push the cloud upward, but the self-gravity 
of the cloud always pulls it down. These up and down 
motions create the complicated structure in the cloud. On 
the other hand, the temperature shows a smooth structure. 
This is due to isothcrmality for each Lagrangian fluid el- 
ement. Only the position of the temperature transition 
region changes in time. 



Figure 3 shows the time evolution of the density. The 
density plots at various times are stacked with time in- 
creasing upward in uniform increments of 0.2to- Be- 
cause the driving force increases linearly with time up 
to t = 10£o > the density changes gradually at first. Af- 
ter t = 10to, the density structure shows many shock 
waves propagating in the cloud, and significant upward 
and downward motions of the outer portion of the cloud, 
including the transition region. After terminating the driv- 
ing force at t = 40to , the shock waves are dissipated in the 
cloud and the transition region falls to around the initial 
position, although it is still oscillating. 



density 




5 10 15 20 25 

Fig. 3. Time evolution of the density. The density plots 
at various times are stacked with time increasing upward in 
uniform increments of 0.2to. 

Figure 4a shows the column density as a function of 
the time averaged position of several Lagrangian fluid ele- 
ments, which are equally spaced at time t = with spacing 



7 



Az = O.lifo starting at z = O.OIHq- The time average is 
calculated between t = 10to and t = 40io, while the driv- 
ing force is input with constant amplitude. The dashed 
line shows the initial distribution. The Lagrangian fluid 
elements have constant enclosed column density as a func- 
tion of time. By using this property, we evaluate the loca- 
tion of Lagrangian fluid elements from the surface density. 
The difference between the initial distribution and that of 
the time average shows that the cloud is lifted up. Figure 
4b shows the time average of the density for each element 
as a function of the time averaged position of each element. 
These time averages were also calculated between t = 10to 
and t = 40to- In contrast to the snapshot in Figure 2, the 
time averaged density structure shows a smooth distribu- 
tion. The dashed line shows the initial density distribution 
of the cloud. The scale height of the time averaged distri- 
bution is about 3 times larger than the initial value. 
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Fig. 4. (a) Column density as a function of the time aver- 
aged position of each Lagrangian fluid element. The dashed 
line shows the initial distribution. The time average is calcu- 
lated between t = 10io and t = 40t . The Lagrangian fluid 
elements have constant enclosed column density as a function 
of time. By using this property, we evaluate the location of 
Lagrangian fluid elements from the surface density. Each La- 
grangian fluid element is equally spaced at time t = with 
spacing Az — 0.1-ffo starting at z = O.Olffo- (b) Time average 
of the density for each fluid element as a function of the time 
averaged position of the element. The dashed line shows the 
initial distribution. 



3.1.2. Velocities in the Cloud 

Figure 5a shows the time evolution of v y at z = 0. Ini- 
tially, it is oscillating about zero sinusoidally with the fre- 
quency of the driving force. However, as time goes on, the 



oscillation shows sharp structures, and it becomes nonsym- 
metric around the mean. The sharp structures are caused 
by nonlinear effects in the cloud such as shock waves. 
Moreover, the mean value shows a deviation from zero 
as time goes on. This implies that the net y-momentum 
of the cloud is not zero and the cloud has a mean motion 
in the y-direction as well as an oscillatory motion. 
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Fig. 5. Transverse velocities versus time. (a) The y- 
component of the velocity at z = as a function of time, (b) 
The mean y-component of velocity of the cloud as a function of 
time, (c) The time average of the mean {/-component of veloc- 
ity for each cycle of the sinusoidal period of the driving force 
as a function of time, (d) The y-component of the velocity 
at 2 = minus the time average of the mean {/-component of 
velocity, as a function of time. 

This mean motion is ultimately caused by the driving 
force. Prior to t = I0£o, the driving force is not symmetric 
in time, because it is proportional to t sin(27Woi). How- 
ever, even if the driving force is symmetric in time, the 
net momentum in the cloud never remains exactly zero 
if we use a driving force like the one in equation (31) 
which is confined to a region near the midplane. The 
net momentum of the cloud is generated by the nonlinear 
and nonsymmetric propagation of this disturbance into the 
stratified cloud, the nonsymmetric restoring forces in the 
system, and nonsymmetric reflection at the cloud bound- 
ary. As disturbances propagate, the restoring force due 
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to the magnetic field is not exactly symmetric in space 
and time. Also, some of the wave momentum is reflected 
at the boundary between the cold cloud and hot external 
medium, and part of it escapes from the cloud. Hence, 
a net y-component of momentum can and does appear 
within the cloud, creating a mean motion. We note that if 
we were modeling a two-dimensional cloud, i.e., the system 
had finite extent in the y-direction, we could reduce the 
net drift with a y-distribution of driving force such that the 
net driving force is zero in the y-direction. Also, if the sys- 
tem was closed in the z-direction, as in a periodic bound- 
ary system, it is possible to add momentum in a regulated 
way so that the net transverse momentum remains exactly 
zero (e.g., Gammie & Ostriker 1996). However, the com- 
bined effects of our 1.5-dimcnsional approximation, cloud 
stratification, and open cloud boundary make it difficult 
to regulate the net drift of the cloud in the y-direction. 

Therefore, for our analysis, we divide v y into two parts. 
The first is the mean velocity which shows the mean mo- 
tion of the entire cloud. The second is the oscillating 
component of the velocity. We calculate the mean y- 
component of velocity of the cloud as 



Io f(t) P v v dz 



»(*) 



(41) 



In this equation, Zf(t) is the full mass position of the cloud, 
which is defined by 



I 



*/(*) E 
pdz = 0.998-^, 



(42) 



and Zf(t) corresponds to the position of the Lagrangian 
fluid element which is initially located at ~ z c , the initial 
position of the transition region of the temperature. The 
time evolution of v m is shown in Figure 5b. 

Figure 5b shows that the mean velocity is still oscillat- 
ing while the driving force is input (t = — 40to)- In order 
to remove the oscillation, we take time average of v m for 
each cycle of the sinusoidal period (to) of the driving force 
and define (v m ) in the following manner. For example, 
(v m ) between nto and (n + l)t is calculated as 

I r-(n+l)t Q 

(v m )(t = nt -(n + l)t ) = - / v mcan (t')dt' , (43) 

l Jnto 

where n is an integer. This calculation is done from n = 
to n — 59. Figure 5c shows ( function of time. Ac- 

cording to the definition, (v m ) has the same value between 
nto and (n + l)to- 

Finally we define the oscillating component of the y- 
velocity as 



= v y - (v m ). 



(44) 



Figure 5d shows the time evolution of v' y . In contrast to 
Figure 5a, it is oscillating around v' y = 0. 

Figure 6 shows the time average of various velocities for 
each Lagrangian fluid element plotted versus the time av- 
erage of the position of each element. Open circles show 
the time averaged root mean square of v Zl filled circles 
show that of v' y , and the dashed line shows the time aver- 
age of the sound speed. The solid line shows the root mean 



square of the time average of the effective one-dimensional 
velocity dispersion 



(45) 



for each fluid element. The time average is taken between 
t = 10t and t = 40£ - This figure shows that the y- 
component of the velocity is the dominant component in 
the cloud. The averaged ^-component of the velocity in- 
creases as a function of (z)t, except for near the midplane. 
This means that the largest velocity dispersion occurs in 
the low density region. This is a similar tendency to that 
of linear Alfven waves. If we assume the WKB approxima- 
tion and no wave dissipation, the energy flux of the waves 
is constant, i.e., 



p (v y) 2 Va = constant, 



(46) 



where Va is the Alfven velocity of the background mag- 
netic field, i.e., 

Va - tt^- (47) 



(4 7 rp) 1 /2- 



This leads to 



-1/4 



(48) 

Although this relation is not exactly applicable to our non- 
linear result, our result shows a similar tendency. We dis- 
cuss this further in the next section. 
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Fig. 6. Time averaged root mean square of various veloci- 
ties for each Lagrangian fluid element versus the time average 
of the position of each element. Open circles show the time 
averaged root mean square of v zi filled circles show that of v' y , 
and the dashed line shows the time average of the sound speed. 
The solid line shows the root mean square of the time average 
of the effective one-dimensional velocity dispersion. 

3.1.3. Pressures in the Cloud 

Figure 7a shows the time averages of thermal pressure, 
magnetic pressure, and dynamic pressures for each La- 
grangian fluid element as a function of the time average of 
the position of each element. The thick dashed line shows 
the thermal pressure, the thick solid line shows the mag- 
netic pressure of the y-component of the magnetic field, 

B 2 

i.e., g^r, the dash-dotted line shows the dynamic pressure 
for the y-component of velocity, i.e., ^p(v' y ) 2 , and the dot- 
ted line shows the dynamic pressure of the z-component 
of velocity, i.e., hpv 2 . The thin straight line shows the 
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magnetic pressure of background magnetic field, and 
the thin dashed line shows the initial thermal pressure. 
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Fig. 7. (a) Time averages of thermal pressure, magnetic pres- 
sure, and dynamic pressures for each Lagrangian fluid element 
as a function of the time averaged position of each element. 
The thick dashed line shows the thermal pressure, the thick 
solid line shows the magnetic pressure of the {/-component of 
the magnetic field, the dash-dotted line shows the dynamic 
pressure of the {/-component of velocity, and the dotted line 
shows the dynamic pressure of the z-component of velocity. 
The thin solid line shows the magnetic pressure of background 
magnetic field, and the thin dashed line shows the initial ther- 
mal pressure, (b) The dynamic pressure for the {/-component 
of velocity (dash-dotted line) and the magnetic pressure of the 
{/-component of the magnetic field (solid line) as a function of 
density. The dotted line shows the scaling expected from the 
WKB theory. 

In contrast to the WKB theory, the dynamic pressure 
of the z-component of the velocity is nonzero in our re- 
sults, although it is smaller than the dynamic pressure 
of the y-componcnt of the velocity. According to the 
small-amplitude theory, the dynamic pressure of the y- 
component of the velocity and the magnetic pressure obey 
cquipartition, i.e., 



1 R 2 



(49) 



Figure 7a shows that such an equipartition is almost sat- 
isfied between z = IHq and z = 5Hq, although there are 
deviations near the midplane. This deviation comes from 
the effect of a nonzero driving force at z — and the sym- 
metrical boundary condition B y (z = 0) = 0. For z > 5H 0} 



a region which contains a small fraction of the total cloud 
mass (see Figure 4), the two energies have distinct spa- 
tial profiles. Near the interface between the cold and hot 
material (the cloud boundary), the magnetic pressure de- 
creases more rapidly than expected from the WKB theory. 
We believe this is due primarily to a standing wave that is 
set up in the outer cloud, with the cloud boundary acting 
as a node for B y and an antinode for v' y . Although many 
different wave modes are generated by the turbulence in 
the cloud, only those which satisfy the boundary condi- 
tion for a standing wave will interfere constructively upon 
reflection. A transverse standing wave can be set up even 
though the boundary itself is moving, since the Alfven 
speed is much greater in the outer cloud than the verti- 
cal speed of the boundary. The influence of the boundary 
reaches well inside the cloud, since the effective wavelength 
of Alfven waves (scaling as p~ x l 2 for waves of fixed fre- 
quency) is also quite large in the low density region. For 
example, at z = 5Ho, where p ~ O.lpoj Alfven waves of the 
input frequency vq = c s o/Hq have wavelength A ~ 4.5Hq; 
at z = 10H , it increases to A ~ 8H - An important result 
is that even though the waves are quasilincar in the outer 
cloud (By <C Bo), the transverse velocity amplitude is 
much larger than expected from equipartition arguments. 

In Figure 7b, we show the dynamic pressure of the y- 
component of velocity and the magnetic pressure of the 
{/-component of the magnetic field as a function of den- 
sity in order to compare with the WKB model predictions. 
We find that the energy in transverse velocities, 1 /2p(v' y ) 2 , 

scales approximately as p 1 / 2 in the outer cloud, the same 
as the WKB theory, although there is a noticeable upward 
turn at z ~ 9Hq due to the standing wave effect. However, 
the wave magnetic energy B 2 /(8ir), which is directly re- 
sponsible for vertical support, decreases more rapidly than 
in the WKB model in the outer cloud, scaling approxi- 
mately as p. We note that strict adherence to the WKB 
predictions is not expected in our model due to nonlinear- 
ity, wave dissipation, and the low frequency v a = c s0 /H of 
the input turbulence. 

3.1.4. Energies in the Cloud 

Figure 8 shows the time evolution of various energies in 
the cloud. These are calculated as follows: kinetic energy 
of the z-component of velocity, 



/• z /(*) i 
E kz {t) = J -^pv 2 z dz; 

kinetic energy of the {/-component of velocity, 

rz/(t) I 



E ky {t) = J -^P V l dz - E krn, 



(50) 



(51) 



where E km is the kinetic energy of the mean motion of the 
cloud, i.e., 



Ekrn (t) 



:/(*) 



pdz; 



(52) 



magnetic energy of the {/-component of magnetic field, 



/■*/« B 2 
E m (t) - / gjrfz; 



(53) 
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sum of the above terms, 



E T (t) = E kz (t) + Ekyit) + E m (t). 



(54) 



In the above equations, the integration was done from 
to Zf(t) because we are interested in energies of the cold 
material. Strictly speaking, the total energy in each case is 
twice the value we calculate due to the mirror symmetric 
boundary condition at z = 0. Equation (51) shows that 
the mean kinetic energy of the cloud is subtracted from 
the kinetic energy of y-component of velocity, so that Ek y 
means the kinetic energy of the oscillating velocity. 
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Fig. 8. Time evolution of various energies in the cloud. The 
thick solid line shows Ek y , the dotted line shows Ekz, the thin 
solid line shows E m , and dash-dotted line shows Et- The val- 
ues of each energy are smoothed out over every one cycle of the 
driving force in order to remove periodic oscillations originating 
from the driving force. 

In Figure 8, the thick solid line shows E^y, the dot- 
ted line shows Ek z , the thin solid line shows E m , and 
dashed-dotted line shows Et . The values of each energy 
are smoothed out over every one cycle of the driving force 
in order to remove periodic oscillations originating from 
the driving force. Among the energies, E^ y and E m axe 
comparable to each other, but Ekz is significantly smaller 
than the others. The sum of energies Et is almost con- 
stant in a logarithmic scale while the input energy is con- 
stant (t = 10<o — 40io), though it has some fluctuations. 
After the driving force is terminated at t = 40io, the en- 
ergies decrease almost exponentially. The decreasing time 
is ~ 8.5io. 

3.2. Parameter Dependence on the Strength of the 
Driving Force 

3.2.1. Density, Velocity and Pressure 

In this paper, we study the effect of changing the 
strength of the driving force, by changing in equation 
(31). 

Figure 9 shows the time evolution of densities for hd = 
20 and 5,^ = 40. The y- velocities at z — 0, v' y , are also 
shown in the figure. This figure shows that a stronger 
driving force causes a larger turbulent velocity, which re- 
sults in a more dynamic evolution of the molecular cloud, 
including stronger shock waves, and larger excursions of 



the cloud boundary. However, after terminating the driv- 
ing force at t = 40£o, the shock waves dissipate and the 
clouds shrink in both cases. 
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Fig. 9. Evolution of clouds with a,d = 20 and dd = 40. (a) 
Time evolution of densities for hd = 20. (b) The oscillating part 
of y- velocity at z = 0, function of time for ad = 20. 

(c) Time evolution of densities for hd = 40. (d) The oscillat- 
ing part of {/-velocity at z = 0, function of time for 
hd = 40. The density plots at various times are stacked with 
time increasing upward in uniform increments of 0.2to- 



Figure 10 shows the time average of the density and ve- 
locities for each Lagrangian fluid clement for hd — 20 and 
hd = 40. The time average is taken between t = 10io 
and t = 40to. This figure also shows that a stronger driv- 
ing force causes the cloud to move further outward and a 
larger velocity dispersion within the cloud. 
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Fig. 10. Density and velocity dispersions for the cases fid = 20 
and &d — 40. (a) Time average of the density for each element 
as a function of the time averaged position of each element for 
fid — 20. The time average is calculated between t = 10io and 
t — 40io- Each Lagrangian fluid element is equally spaced at 
time t = with spacing Az — O.IHq starting at z = 0.01-ffo- 
(b) The time averaged root mean square of various velocities for 
each Lagrangian fluid element plotted versus the time average 
of the position of each element for a d = 20. Open circles show 
the time averaged root mean square of v z , filled circles show 
that of v' y , and the dashed line shows the time average of the 
sound speed c s . The solid line shows the root mean square of 
the time average of the effective one-dimensional velocity dis- 
persion, (c) The time average of the density for each element 
as a function of the time averaged position of each element for 
fid = 40. (d) The time averaged root mean square of various 
velocities for each Lagrangian fluid element plotted versus the 
time average of the position of each element for fid = 40. 

Figure 1 1 shows the time average of pressures for dd = 
20 and da — 40. When the driving force is strong, the mag- 
netic pressure and dynamic pressure of the z-componcnt 
of velocity become significantly larger than the thermal 
pressure. 
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Fig. 11. Time averages of the thermal pressure, magnetic 
pressure, and dynamic pressures for each Lagrangian fluid el- 
ement as a function of the time averaged position of each el- 
ement. The upper panel is for fid = 20 and the lower panel 
is for fid = 40. The thick dashed line shows the thermal pres- 
sure, the thick solid line shows the magnetic pressure of the 
y-component of the magnetic field, the dash-dotted line shows 
the dynamic pressure of the y-component of velocity, and the 



dotted line shows the dynamic pressure of the ^-component of 
velocity. In both panels, the thin solid line shows the magnetic 
pressure of background magnetic field, and the thin dashed line 
shows the initial thermal pressure. 

3.2.2. Energy 

Figure 12 shows the time evolution of Et for both 
dd — 20 and a d = 40. These values are also smoothed 
out over every one cycle of the driving force to remove 
periodic oscillations. In the case of dd = 20, the energy 
is almost constant in a logarithmic scale while the driv- 
ing amplitude is constant (t = 10t — 40to). However, 
in the case of dd = 40, the energy is still gradually in- 
creasing until t = 30io, but afterwards becomes almost 
constant until t = 40to- After terminating the driving 
force at t = 40to, both energies decrease almost exponen- 
tially. The energy decreasing times, td, for each parameter, 
which are estimated by fitting an exponential function, are 
listed in Table 1. In this study, we terminate the driving 
force at t = 40io in every case for simplicity. However, 
we found that the energy decreasing time can vary some- 
what depending on when we terminate the driving force. 
Accounting for this as well as the fitting error of the expo- 
nential function, we conclude that the energy decreasing 
times have a range of variation about ±2i of the values 
listed in Table 1. 
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Fig. 12. Time evolution of Et for both fid = 20 and fid = 40. 
The values of each energy are smoothed out over every one cy- 
cle of the driving force in order to remove periodic oscillations 
originating from the driving force. 

3.2.3. Correlations between Velocities and Sizes 

Here we investigate the correlation between the velocity 
dispersion and the height of the cloud. Figure 13a shows 

the time averaged velocity dispersions {a 2 ) X J 2 of different 
Lagrangian fluid elements for different dd, as a function 
of (z) t - The open circles correspond to Lagrangian fluid 
elements whose initial positions are located at z = 2.51, 
which is close to the edge of the cold cloud. The filled cir- 
cles corresponds to Lagrangian fluid elements whose initial 
positions are located at z — 0.61, which is approximately 
the half-mass position of the cold cloud. Each circle cor- 
responds to a different value of dd- These values are sum- 
marized in Table 1. The dotted line shows 



0.5 
t ■ 



(55) 
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This figure shows that the velocity dispersions have a good 
correlation with the heights of the molecular clouds. 
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Fig. 13. Global properties of an ensemble of clouds with 
different driving strengths 2^. (a) Time averaged velocity dis- 
persions of different Lagrangian fluid elements for different 3d, 
as a function of time averaged positions. The open circles cor- 
respond to Lagrangian fluid elements whose initial positions 
are located at z — 2.51, which is close to the edge of the cold 
cloud. The filled circles correspond to Lagrangian fluid ele- 
ments whose initial positions are located at z = 0.61, which is 
approximately the half-mass position of the cold cloud. The 
dotted line shows {a 2 ) x J 2 cx (z)t' 5 ■ Each circle can be associ- 
ated with a particular model in our study by comparison with 
the numbers in Table 1. (b) Time averaged velocity disper- 
sions as a function of the mean Alfven velocity of the cloud. 
The dotted line shows (o 2 ) 1 ^ 2 oc Va- 

Figure 13b shows the correlation between the velocity 
dispersion and mean Alfven velocity of the cloud Va at 
the mean position (z) t of the Lagrangian elements. The 
dotted line shows 

{v 2 ))' 2 cx V A . (56) 

This figure shows that the velocity dispersions have a good 
correlation with the mean Alfven velocity, defined by 



where 



V A 



P = 



B 



Tip 



(57) 



(58) 



is the mean density and £ is the column density for each 
Lagrangian element. We discuss the meaning of these cor- 
relations in the next section. 



4.1. 



4. DISCUSSION 

Velocity Dispersion and Equilibrium 

In Figure 13a, we found that the velocity dispersion 
obeys the lincwidth-size relation 



(59) 



This relation is satisfied (at both the full-mass and half- 
mass Lagrangian positions) for an ensemble of the clouds 
with different strengths of the driving force. It is consis- 
tent with observational results of molecular clouds (Lar- 
son 1981; Myers 1983; Solomon et al. 1987). Within any 

individual cloud, (cr 2 )]^ 2 also increases toward the cloud 
boundary. 

A similar but not identical relation, relating the veloc- 
ity dispersion at z = to the size of the cloud, is ex- 
pected in virtually any one-dimensional model, regardless 
of the form of spatial variation of the pressure within the 
cloud. This is because an integral over the vertical force- 
balance equation reveals that the total pressure at z = 0, 
-£tot,o must equal the weight of the accumulated gas above, 
7rGE 2 /2, assuming that the surface pressure is negligible. 
If ftot,o = Po< j2 s.o7 where a e fi,o is the effective velocity dis- 
persion at z = 0, and £ = 2poZ, in which Z is the typical 
size scale of the cloud, then one naturally obtains 



O~eS,0 oc 



z l/2 



(60) 



for an ensemble of clouds of fixed total column density 
£ but varying cr e ff,o- This relation applies to the Spitzer 
equilibrium state, in which <7 e ff,o = c s and Z — H n . It also 
applies to the equilibrium state calculated by Martin ct al. 
(1997), in which cT e ff,o is the effective velocity dispersion 
associated with Alfven waves at z = 0. While a similar 
relation will inevitably apply to our nonlinear model as 
well, we note that the relation (59) is more general, in 
that it relates the velocity dispersion at the position (not 
z = 0) of a Lagrangian mass element to the position of 
that element. This more general correlation is a less pre- 
dictable property of our time-averaged equilibrium state. 
It is also closer to what is often measured, since optical 
depth effects may mean that a measured velocity disper- 
sion samples the largest scales of an observed cloud rather 
than the center. 

Equation (59) is strongly related to another relation 



{o 2 )\ 12 oc V A 



(61) 



which is shown in Figure 13b, once more for both the half- 
mass and full-mass positions for our ensemble of clouds. 
This relation is also consistent with observational results of 
molecular clouds (Crutcher 1999; Basu 2000). For exam- 
ple, Figure lb of Basu (2000) shows an excellent correla- 
tion between the line-of-sight component of the large-scale 
magnetic field B\ a& and ap 1 ^ 2 for observed clouds of widely 
varying values of p, Bi os , and a, essentially the same cor- 
relation as equation (61) if Bi os cx B on average. 
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Here we would like to point out that the relation (61) 
does not necessarily imply that turbulent motions are due 
to Alfven waves, although that is primarily the case in 
our simulation. The relation is far more general in that 
it is obtained when the velocity dispersion is caused by 
any mechanism which results in the clouds being in virial 
balance between gravitational and turbulent energies, as 
well as having a mass-to-flux ratio that is approximately 
uniform from cloud to cloud (Mouschovias 1987; Shu et al. 
1987). Rather than imagining a time- independent mean 
Alfven speed in a cloud (strictly true in an incompress- 
ible medium or a periodic box simulation with a fixed 
mean density), and turbulent motions becoming compa- 
rable to that speed, our simulation reveals that it is the 
mean Alfven speed of any cloud which readjusts to a new 
value as p drops to accommodate a particular level of tur- 
bulent driving. In this view, the turbulent dispersion a is 
the more fundamental quantity, dependent on the particu- 
lar source of driving, and Va is a quantity which readjusts 
to become comparable to a. 

Finally, we note that the relation (61) depends on the 
definition of the mean Alfven velocity. If we use the local 
Alfven velocity of each particle, we could not get a clear 
correlation like equation (61). The mean Alfven velocity 
is often calculated observationally using the mean density, 
as in equations (57) and (58), because it is difficult to di- 
rectly measure a local Alfven velocity. The local Alfven 
velocity in a gravitationally stratified cloud often takes on 
values much different than the overall mean quantity Va 
used in equation (61). For example, the time averaged lo- 
cal Alfven velocity near the edge of the cloud is about 2.6 
times greater than Va for the case of ad = 30. 

4.2. Relation to Linear Model and Chandrasekhar- Fermi 

Formula 

In our simulation, the velocity component parallel to the 
background magnetic field (v z ) is generated by the nonlin- 
ear effect of the waves. However, the time-averaged magni- 
tude of v z is significantly less than the time-averaged mag- 
nitude of v' y in all our models (see Figures 6, 7, 10, and 11). 
Although the waves are nonlinear, the coupling between 
the Alfven and slow mode MHD waves is not so strong as 
to destroy an approximate equipartition between By/ (Sir) 
and 1 /2p{v' y )' 2 throughout much of the cloud, which is the 
expected result for purely transverse linear Alfven waves. 
However, this equipartition does break down in the outer 
part of the cloud, where standing wave motions are estab- 
lished in which B y has a node at the cloud boundary, and 
v'y has an antinode. 

The breakdown of the v' y v.s. B y relation in outer 
cloud means that the original Chandrasekhar-Fermi for- 
mula breaks down there as well. If we assume that the 
dispersion of polarization angle of the magnetic field (58) 
is related to the velocity dispersion as 



Bo Va 



(62) 



we can estimate the strength of the magnetic field of the 
cloud by using the observational values of polarization an- 
gle, density, and velocity dispersion. This yields 



where a is a nondimcnsional factor and equals 1 for linear 
Alfven waves. The relation (63) with a = 1 was proposed 
by Chandrasekhar & Fermi (1953) in order to estimate the 
strength of the background magnetic field in the interstel- 
lar medium. The generalized form with a =/= 1 can be fit 
to our simulation if we define a in our simulation as the 
square root of the ratio of time averaged magnetic pressure 
to dynamic pressure for each Lagrangian fluid element, i.e., 



(pK) 2 /2) t 



1/2 



(64) 



The resulting distribution of a as a function of (z) t for 
the case of a<j = 30 is shown in Figure 14. The time 
averaged density, which was shown in Figure 4, is also 
shown as a dashed line to clarify the edge of the cloud. 
It shows that a is close to unity inside the cloud in the 
region where most of the mass is enclosed (see Figure 4), 
but it decreases near the surface of the cold cloud, where 
the standing wave effect becomes important; the minimum 
value is about a — 0.23. Outside of the cloud {{z) t > 12H 
in Figure 14), a — 1 since the waves are linear there due 
to the low density and high ambient Alfven speed. There- 
fore, while our nonlinear model supports the use of the 
Chandrasekhar-Fermi formula throughout most of a strat- 
ified cloud, we caution against its use with a = 1 near the 
surface of a cloud. 
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Fig. 14. Parameter a which appears in the Chandrasekhar- 
Fermi formula (eq. [63]), as a function of the time averaged po- 
sition of fluid elements within the cloud for the case of hd = 30. 
The time averaged density, which was shown in Figure 4, is 
also shown as a dashed line, mainly to clarify the position of 
the edge of the cloud. 

4.3. Global Oscillations 

Figure 15 shows the time evolution of the position of a 
Lagrangian fluid element whose initial position is z = 2.51, 
for the standard model with driving strength aa = 30. 
This corresponds to the motion of the outer edge of the 
cloud. The motions resemble longitudinal normal mode 
oscillations, with the excursions of the outer cloud very 
similar to free-fall. The dotted line in the figure shows 
the trajectory for free-fall motion of the fluid element for 
several different time intervals. The peak of the trajecto- 
ries are chosen to coincide with a peak of the oscillation in 
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the computed model. The paths are parabolic due to the 
constant gravity acting on a comoving mass shell in the 
one-dimensional approximation. Similar to the transverse 
standing wave that is set up in the outer cloud, the lon- 
gitudinal motions in this region also resemble a standing 
wave pattern with an antinodc of v z at the cloud bound- 
ary. The outermost part of the cloud suffers the greatest 
displacements, as in the case of a pulsating star. Once 
the internal driving is discontinued (t > 40 to), the outer 
surface also moves inward in nearly a free-fall manner. 

20 f ' ' ' 1 ' ' ' 1 ' ' ' ] 




20 40 60 

t/to 

Fig. 15. Time evolution of the position of a Lagrangian fluid 
element whose initial position is located at z = 2.51 for the case 
of &d = 30. The dotted line shows the trajectories of free-fall 
motion for the element. 

Figure 3 and Figure 15 also show that some residual os- 
cillations remain on the largest scale even after the damp- 
ing of most of the internal turbulent energy. This effect is 
reminiscent of the recent observation of Lada et al. (2003) 
that implied a global oscillation of Barnard 68, a cloud 
with only subsonic internal turbulence. We hypothesize 
that Barnard 68 may have dissipated most of the internal 
turbulence left over from its formation, but that only a 
large-scale oscillation remains, as in the late stage of our 
simulation after turbulent driving is discontinued. How- 
ever, we caution that ours is a one-dimensional model, 
and globally coherent standing-wave motions need to be 
demonstrated in a multi-dimensional model. 

4.4. Dissipation of Energy 

The dissipation of energy that we input into nonlin- 
ear transverse Alfven modes is caused though either shock 
waves or grid-scale dissipation, which is similar to the re- 
sults of Stone et al. (1998). In addition to the dissipa- 
tion, a part of the energy escapes from the cloud in our 
simulation. We measured the Poynting flux, which is the 
dominant energy flux, at the full-mass position and the 20- 
percent mass position, which is just outside the region we 
input the driving force. The time integral of the Poynting 
flux at the full-mass position is about 30 percent of that 
at the 20-percent mass position. 

In our model, the majority of dissipation occurs due to 
the transfer of energy to small scales via nonlinear steep- 
ening and/or a turbulent cascade, followed by damping on 
the grid scale though numerical resistivity. This is because 
the kinetic energy of longitudinal motions, which are the 
cause of shock waves, is smaller than that of transverse mo- 
tions. The primary dissipation of Alfvcnic motions then, is 
not due to coupling to slow modes, although this is present 



at a significant level. This is similar to the finding of the 
one-dimensional numerical model of Gammie & Ostrikcr 
(1996), although the same authors claim that in a multi- 
dimensional simulation, compressive effects become a ma- 
jor source of dissipation (Stone et al. 1998); however, see 
Cho & Lazarian (2003), who find that the primary dissi- 
pation mechanism in a three-dimensional simulation is not 
the coupling of Alfven modes to slow and fast modes. 

So if the primary source of dissipation in this simula- 
tion is grid-scale dissipation after nonlinear steepening of 
waves (and perhaps some cascade to smaller wavelengths), 
then is the dissipation rate a physical effect or an artifact 
of the simulations? To answer this, we have increased the 
resolution of our numerical simulations until the dissipa- 
tion rate is nearly independent of grid size for the case of 
the typical parameter. Hence, the numerical diffusion is 
no longer resolution-determined at our current resolution 
of 50 points per initial scale height in the region where 
most of the mass is located. If we use a smaller grid size, 
the final scale of the dissipation is smaller, but the rate of 
dissipation stays nearly the same. In reality, there must be 
a physical effect such as ion-neutral friction which would 
damp the waves at small scales. 

The dissipation time of our results are a few crossing 
times of the time-averaged scales of the clouds (see Table 
1). They are a bit longer than estimated from periodic box 
simulations. We think that this comes from the generation 
of longer-wavelength modes as the waves travel to low den- 
sity regions near the cloud's surface. However, it is already 
known that one-dimensional simulations have lower dissi- 
pation rates than two or three-dimensional ones (Ostriker 
et al. 2001). Therefore, higher dimensional global simula- 
tions will give the final answer in the future. 

4.5. Future Work 

This is the first in a series of papers. In the next pa- 
per, we will conduct a complete survey of the effect of im- 
portant parameters such as (3q and i>q. Additionally, the 
study of random, rather than sinusoidal, disturbances will 
be considered. Moreover, we will study the case of circu- 
larly polarized Alfven waves by including an x-component 
of motions in the simulation. A circularly polarized wave 
is possibly less dissipative than a linearly polarized wave, 
since in the theoretical limit of an infinite wave train, a 
flux of circularly polarized waves has no associated mag- 
netic pressure gradient and thereby induces no compres- 
sion of gas. Although conversion of Alfven modes to slow 
modes is already not a dominant source of dissipation in 
our model, it will be useful to see if the dissipation rate 
is even less than measured here when circularly polarized 
waves propagate in a finite-sized cloud. 

We also expect to include ion-neutral friction in a fu- 
ture simulation. This effect is expected to be important 
in molecular clouds, especially in this model as a damp- 
ing mechanism of Alfven waves (Kulsrud & Pearce 1969; 
Zweibel & Josafatsson 1983). It would also change the 
density structure of a cloud lifted up by wave pressure, as 
demonstrated in the case of linear Alfven waves by Martin 
et al. (1997). 

Finally, we believe that a multi-dimensional turbulent 
simulation in a gravitationally stratified medium is very 
necessary for the future, especially to obtain more defini- 
tive dissipation rates in a stratified cloud and to study 
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the feedback to turbulence due to gravitational contrac- 
tion and differential rotation. 

5. SUMMARY 

We have performed a numerical simulation of nonlin- 
ear MHD waves in a stratified molecular cloud. Our main 
results are as follows: 

1. Due to the effective pressure of MHD turbulence, 
our one-dimensional cloud is lifted upward and es- 
tablishes a steady-state characterized by oscillations 
about a time-averaged equilibrium state. The outer, 
low density parts of the cloud undergo the largest os- 
cillations, which have the character of free-fall mo- 
tions. After turbulent driving is discontinued, the 
cloud falls back toward the initial equilibrium, but 
some large-scale oscillations of the cloud are the 
longest-lived modes. 

2. The nonlinear effect of the wave propagation results 
in a significant conversion of energy from Alfvcn 
waves to compressive slow (acoustic) modes and 
consequent shock formation. However, the energy 
in transverse modes always remains significantly 
greater than that in longitudinal modes. Further dis- 
sipation of transverse modes takes place through the 
generation of smaller scale structure through nonlin- 
ear steepening or a cascade, culminating in resistive 
dissipation on the grid-scale. 

3. Within each cloud, the magnetic pressure associated 
with wave motion, (By/(8n))t, is in approximate 
equipartition with the kinetic energy in transverse 
motions (p(v' y ) 2 /2) t in the region containing most of 
the mass. However, in the low density outer regions, 
the wave amplitudes have the character of stand- 
ing waves, such that B y has a node and v' y has an 
antinode. All of this means that the Chandrasekhar- 
Fermi model applied to a stratified cloud will have 
a multiplicative constant a which is « 1 with weak 



spatial dependence through much of the cloud, but 
which drops to much lower values near the edge of 
the cloud. 

4. After turbulent driving is discontinued, the dissipa- 
tion time of the cloud turbulence is ~ 10 to or several 
crossing times of the time-averaged equilibrium state 
during turbulent driving. These times are up to a few 
times longer than those found in multi-dimensional 
periodic simulations. This would be due partly to 
fewer dissipation avenues in a one-dimensional ap- 
proximation, but also partly to the generation of 
long-wavelength modes due to the cloud stratifica- 
tion. 

5. For an ensemble of clouds with different levels of in- 
ternal driving, we find the relation a oc Z 5 , where a 
is the time averaged one-dimensional velocity disper- 
sion measured at a size scale Z of the time averaged 
equilibrium state of the cloud. 

6. For the same ensemble of clouds, the above relation 
also implies that a oc Va, where Va is the mean 
Alfvcn velocity determined from the mean density 
of the cloud. This relation is in agreement with 
observed properties of magnetized clouds and cloud 
cores. 
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Table 1 

Dissipation time as a function of a d - Time- averaged height and time- averaged velocity dispersion as a 

FUNCTION OF a d FOR TWO LAGRANGIAN ELEMENTS: z(t = 0) = 2.51, APPROXIMATELY THE FULL-MASS POSITION, 
AND z(t = 0) = 0.61, APPROXIMATELY THE HALF-MASS POSITION. 
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